/*
 * Program to generate the initials spacings between the y's
 * Salem Machaka
 * 04/10/98
 */

#include"mle.h"

/*
 * This function Calculates the initial spacings between the 
 * probes before the continuous optimization algorithm starts.
 * This is used as a "good" starting point in order not to
 * spend too much time minimizing.
 * This is solved as a linear programming problem that tries
 * to minimize the initial spacing for the probes that have 
 * a maximal number of clones that hybridize to them
 */

void Initial_Ys(int *order, double *y, int nprobe, int nclone,
        int M, int gap, int **joint) {
    
    int w, mmax, mmin, msum, j, un, t1, t2, *link, *ulink;
    double A, eps, t3, t4, eps1, sum;


    eps1 = 0.1;
    msum = 0;
    mmax = 0;
    mmin = nclone;
    w = 0;

    link = Int_Vector(1, nprobe);
    ulink = Int_Vector(1, nprobe);
    un = 0;

    /*
     * w = number of kinked probes
     * msum = total sum of linking clones for all probes
     * mmax = Maximum # of linking clones between any 2 probes
     * mmin = Minimum # of linking clones between any 2 probes
     */
    for (j = 2; j <= nprobe; j++) {
        t1 = joint[order[j]][order[j - 1]];
        if (t1 > 0) {
            w++;
            msum += t1;
            if (mmax < t1) mmax = t1;
            if (mmin > t1) mmin = t1;
            link[w] = j;
        } else {
            un += 1;
            ulink[un] = j;
        }
    }

    /*
     * A = space for all linked probes
     */
    A = (double) gap - ((nprobe + 1 - w)*(M - eps1));
    /*
     * If A<=0, no more space is left for the linked
     * probes. Unlinked probes = gap/(nprobe-w+1).
     * Linked probes = 0
     */
    if (A <= 0.0) {
        y[1] = gap / (double) (nprobe + 1 - w);
        for (j = 1; j <= un; j++)
            y[ulink[j]] = y[1];
        for (j = 1; j <= w; j++)
            y[link[j]] = 0.0;
    }        /*
   * Else , run linear programming model
   */
    else {
        eps = A / (2 * w);
        y[1] = (gap - A) / (nprobe + 1 - w);
        for (j = 1; j <= un; j++)
            y[ulink[j]] = y[1];
        if (((M - eps) * w) <= A) {
            for (j = 1; j <= w; j++)
                y[link[j]] = M - eps;
        } else {
            t1 = (mmax * w) - msum;
            t2 = msum - (mmin * w);
            if (t1 == 0) /* which implies that t2 would be zero as well */ {
                for (j = 1; j <= w; j++)
                    y[link[j]] = A / w;
            } else {
                t3 = (A - (eps * w)) / t1;
                t4 = ((w * (M - eps)) - A) / t2;
                if (t3 < t4) {
                    for (j = 1; j <= w; j++) {
                        un = joint[order[link[j]]][order[link[j] - 1]];
                        y[link[j]] = ((A * (mmax - un))+(eps * ((w * un) - msum))) / t1;
                    }
                } else {
                    for (j = 1; j <= w; j++) {
                        un = joint[order[link[j]]][order[link[j] - 1]];
                        y[link[j]] = (((M - eps)*(msum - (w * un)))+(A * (un - mmin))) / t2;
                    }
                }
            }
        }
    }

    /*
     * Assigning initial value to y[n+1]
     */
    sum = 0.0;
    for (j = 1; j <= nprobe; j++)
        sum += y[j];
    y[nprobe + 1] = gap - sum;

    Free_Int(link, 1);
    Free_Int(ulink, 1);


}

/*
 * Function to evaluate the constant vector R
 */
void Get_R(double *R, int N, int nprobe, int nclone, int M, double **aa, int *order) {

    int i, j;
    double sum;

    for (i = 1; i <= nclone; i++) {
        sum = 0.0;
        for (j = 1; j <= nprobe - 1; j++)
            sum += aa[i][order[j]] * aa[i][order[j + 1]];
        R[i] = N - nprobe * M + M*sum;
    }

}

